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Preface 


One  limit  to  the  information  retrievable  from  an  image  Is  the 
"diffraction  limit”  which  translates  to  limited  resolution.  Superresolved 
Imaging  is  a  field  devoted  to  increasing  resolution  beyond  this  limit. 
Whereas  Iterative  superresolution  has  been  shown  to  work  In  general,  it 
traditionally  requires  a  knowledge  of  the  finite  dimensions  of  the  object 
This  thesis  demonstrates  that  the  same  technique  can  be  used  without 
tnis  knowledge.  The  technique  has  also  been  modified  to  use  the  more 
efficient  Hartley  transform  instead  of  the  Fourier  transform. 


As  indirect  contributors  to  the  results  shown  herein,  my  thanks  go 
to  the  following:  my  wife  Sandra,  for  being  her;  Dr.  Theodore  E.  Luke, 
for  the  long  leash;  Dr.  Steven  K.  Rogers  for  his  advice  on  phase 
retrieval;  LtCol  James  P.  Mills  for  challenging  my  work;  Bruce  Hornsby, 
Pink  Floyd  and  Rush  for  inspiration;  and  John  T.  Kelly,  my  father,  for 
always  asking  me  "Why?". 
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ABSTRACT 

/ 

Resolution  In  an  image  can  be  increased  by  an  iterative  technique 
(Introduced  by  Gerchberg)  which  effectively  continues  the  known,  partial 
spectrum  beyond  the  limits  Imposed  by  an  optical  system.  This  increased 
resolution  is  called  superresolution.  Historically,  an  important  con¬ 
straint  assumed  for  this  technique  was  the  knowledge  of  the  finite  dimen¬ 
sions  of  the  object  such  that  the  object  energy  outside  these  dimensions 
must  be  zero.  It  is  shown  in  this  thesis  that  by  a  change  of  object 
space  geometry,  the  semicircular  field  of  view  of  an  optical  system 
provides  a  natural  dimensional  constraint  which  can  be  used  instead  of 
the  object  dimensions  to  achieve  superresolution.  A  further  modification 
of  the  iterative  technique  involves  using  the  Fast  Hartley  Transform 
(FHT)  Instead  of  the  Fast  Fourier  Transform  (FFT).  The  FHT  is  inherently 
faster  and  requires  less  computer  memory  than  the  FFT. 


SUPERRESOLVED  IMAGING  USING  THE  HARTLEY 
TRANSFORM  AND  A  SPHERICAL  OBJECT  SPACE 


I.  INTRODUCTION 


According  to  Gerchberg  (1:711)  and  supported  by  Fienup  (2:161), 
seemingly  lost  spatial  information  about  an  object  might  be 
recoverable  if  constraints  are  known  in  both  the  object  and  object 
spatial  frequency  domains.  The  transform  generally  relating  these 
domains  is  the  Fourier  transform.  One  constraint  in  object  space 
which  is  commonly  used  is  the  finite  dimension  of  the  object. 

However,  if  this  knowledge  Is  gained  from  received  data  (ie  autocor¬ 
relation  data),  then  it  may  be  highly  suspect  in  real  scenarios  where 
detector  noise,  diffracted  light  and  undesireable  objects  in  the  field 
of  view  retard  the  desired  object's  autocorrelation  properties. 

The  recovery  of  lost  spatial  information  is  of  great  importance 
to  both  civilian  and  military  applications  of  imaging.  One  loss 
mechanism  is  derived  from  the  relationship  between  the  extent  of  the 
received  spatial  frequency  function  (the  spectrum)  and  the  maximum 
resolution  obtainable  from  a  given  optical  system  aperture.  As  the 
aperture  size  increases,  higher  resolution  is  possible  due  to  the 
larger  spectrum  passing  through  the  aperture.  However,  there  are 
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many  practical  restrictions  to  aperture  size.  Another  method  of  in¬ 
creasing  resolution  is  to  continue  the  finite  or  "band-limited" 
spectrum  of  a  small  aperture  further  into  frequency  space  to  simulate 
using  a  larger  aperture.  This  is  possible  due  to  the  analytic  nature 
of  an  object's  spectrum  (3:932).  Such  a  continuation  can  be  ac¬ 
complished  by  the  Gerchberg  technique  given  sufficient  constraints. 

The  enhanced  spectrum  can  then  be  transformed  into  an  image  with 
higher  resolution  than  the  image  formed  with  the  initial  spectrum. 

Other  investigations  of  superresolved  imaging  can  be  found  in  appendix 

A. 


The  purpose  of  this  thesis  is  to  examine  an  alternative  con¬ 
straint  in  object  space  and  apply  this  constraint  to  the  iterative  su¬ 
perresolution  technique.  The  new  constraint  naturally  appears  when  a 
spherical  object  space  is  chosen  instead  of  the  traditional  infinite 
rectangular  object  space.  Also,  since  this  technique  is  Iterative,  a 
fast  Hartley  transform  will  be  used  instead  of  the  more  obvious  fast 
Fourier  transform  due  to  the  former's  greater  speed  and  smaller 
memory  requirements.  The  development  will  only  consider  a  one¬ 
dimensional  object  (le  an  object  function  of  only  one  variable). 


II.  BACKGROUND 


Iterative  Superresolution. 


To  paraphrase  Gerchberg  (1:712),  the  superresolution  algorithm  is 
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as  follows.  The  band-limited  spectrum  Is  inverse-Fourier  transformed 
to  create  an  image.  The  energy  which  falls  outside  of  the  known  ob¬ 
ject  dimensions  in  image  space  is  set  to  zero.  This  modified  image 
function  is  Fourier  transformed  to  generate  a  new  spectrum.  The  new 
spectrum  is  replaced  by  the  Initial  band-limited  spectrum  only  within 
that  band-limitation.  Repeating  this  routine  gradually  builds  the  ini¬ 
tial  spectrum  out  beyond  the  limitations  imposed  by  the  optical  sys¬ 
tem  which  generated  it. 

The  explanation  of  the  above  method's  ability  to  continue  the 
finite  spectrum  is  fairly  simple.  The  actual  spectrum  of  any  finite 
object  continues  in  frequency  space  indefinitely.  Thus,  the  band- 
limited  spectrum  can  be  considered  as  the  actual  spectrum  plus  an 
"error"  spectrum  which  is  first,  equal  and  opposite  the  actual 
spectrum  outside  the  frequency  band  passed  by  the  optical  system  and 
second,  zero  inside  this  band  (making  the  total  spectrum  outside  the 
band  zero).  The  linear  properties  of  the  Fourier  transform  allow 
separate  treatment  of  each  constituent  spectrum.  The  actual 
spectrum  obviously  Inverse-transforms  to  the  original  object  function. 
The  error  spectrum,  being  zero  for  a  finite  region  and  generally  non¬ 
zero  elsewhere,  will  inverse-transform  to  an  image  error  function 
which  has  energy  outside  the  object  constraints  in  image  space. 

Since  this  is  the  energy  which  is  zeroed  in  the  algorithm,  then  by 
Parseval's  theorem,  only  the  energy  of  the  error  spectrum  is 
decreased.  As  the  error  spectrum  decreases,  the  actual  spectrum  is 
revealed. 
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Use  of  the  Hartley  Transform. 

As  an  important  note,  the  above  superresolution  technique  Is  not 
necessarily  confined  to  a  Fourier  transform  space  relation.  Such  a 
scenario  is  commonly  found  in  optics  however,  and  justifies  its  initial 
establishment.  In  general,  any  transform/inverse-transform  relation 
which  maps  energy  beyond  constraints  in  one  space  and  has  a  con¬ 
straint  in  the  other  space  should  work.  While  the  Hartley  transform 
provides  one  of  these  relations,  it  will  not  be  introduced  as  such.  A 
Fast  Hartley  Transform  (FHT)  will  be  used  in  this  thesis  only  in  its 
capacity  to  make  the  superresolution  algorithm  more  efficient  in  com¬ 
putation  time  and  computer  memory  requirements.  Bracewell  (4:12)  has 
shown  that  the  Hartley  transform,  h(x),  can  be  determined  from  the 
Fourier  transform,  f(x),  by 

h(x)  =  fr*al(x)  -  fiaag(x)  (1) 

Since  the  ultimate  objective  is  a  superresolved  image,  not  necessarily 
a  continued  Fourier  spectrum,  conversion  back  to  the  Fourier  domain 
is  unnecessary.  Thus  the  superresolution  algorithm  will  be  performed 
from  image  space  to  frequency  space  via  the  FHT,  with  the  frequency 
function  called  the  Hartley  spectrum.  The  algorithm  is  summarized  in 
Figure  one.  The  convergence  criterion  is  simply  a  given  number  of 
Iterations. 
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The  constraint  in  frequency  space  is  now  the  known,  band-limited 
Hartley  spectrum  converted  from  the  Fourier  spectrum.  The  array 

"RESULT"  contains  the  discrete  values  of  either  the  Hartley  spectrum  _• 

or  image  functions,  depending  on  the  step  of  the  algorithm. 
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III.  THEORY 


The  relation  between  the  two  spaces  used  in  the  superresoiution 
algorithm  is  defined  empirically  by  the  propagation  of  spherical  waves 
from  the  object  in  accordance  with  the  Fresnel-Kirchhoff  diffraction 
formula  which  assumes  that  the  separation  between  spaces  is  much 
greater  than  the  wavelength  of  radiation  observed  (5:41).  If  the  ob¬ 
ject  and  optical  system  aperture  dimensions  are  also  much  smaller 
than  this  separation,  then  it  is  common  practice  to  approximate  the 
diffraction  formula,  and  thus  the  space  relation,  with  a  Fourier 
transform.  It  is  important  to  emphasize  that  this  relation  only  ap¬ 
plies  to  a  small  maximum  object  dimension.  Even  though  the  rectan¬ 
gular  object  space  historically  used  extends  Indefinitely,  the  results 
of  an  inverse-Fourier  transform  which  map  outside  the  small  region 
can  not  be  assumed  accurate  because  such  an  operation  violates  the 
geometry. 

There  are  two  reasons  for  choosing  a  spherical  object  space 
over  the  rectangular  object  space.  First,  the  propagation  geometry 
from  such  a  curved  space  should  more  closely  approximate  a  Fourier 
transform  over  the  entire  object  space.  Second,  with  a  spherical 
space,  only  a  hemisphere  of  light  can  enter  a  windowed  aperture  plane. 
This  places  a  natural  constraint  on  the  object  space  which  may  be 


used  to  reduced  error  energy. 


Transform  Formulation. 


! 

I 


c 


The  traditional  diffraction  formula  can  be  written  as 


D(x) 


r/2 

c  j  0(a) 

-  *V2 


exp(ikr) 

r 


Z(x,a,r)da 


(2) 


where  D(x)  is  a  complex  distribution  at  the  system  aperture  plane,  c 
is  a  complex  constant,  0(a)  is  the  object  function,  k  =  2 ir/\  (X  Is 
the  wavelength)  and  Z(x,a,r)  is  the  obliquity  function  (5:58).  The 
geometry  is  shown  in  Figure  2. 


( 
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The  limits  of  integration  are  representative  of  the  semicircular 
geometry.  Of  course  there  may  be  object  energy  outside  this  semi¬ 
circle,  but  none  of  the  energy  which  enters  the  system  aperture  is 
from  this  outside  regime.  Ve  make  an  initial  but  standard  approxima¬ 
tion  that  if  R  is  very  large,  then  r  ~  R  for  all  amplitude  terms. 

This  also  affects  the  obliquity  function  such  that  it  becomes  only  a 
function  ota.  Thus 

*/2 

D{x)  =  ^  j*0(a)expllkr)Z(a)da!  (3) 

-T/2 

We  now  explore  the  phase  term.  The  distance  r*can  be  written 
as  (xo-x)a+(z<>)2  by  the  Pythagorean  theorem.  In  spherical  coordinates, 
Zo  =  Rcos(a)  and  xo  =  Rsin<«)  .  Thus 


We  rename  the  expression  in  brackets  b.  In  examining  the  mag¬ 
nitude  of  b,  with  the  scenario  that  x<<R,  we  find  two  extreme  cases 
depending  on  the  magnitude  of  the  sine  function.  If  sin(cr)  =  1  or  -1 
then  the  magnitude  of  b  is  given  by  (x/R)(x/R  -  2),  whereas  if  sin(ar)  = 
0  then  this  magnitude  becomes  (x/R)2.  In  both  cases,  the  magnitude  of 


b  is  less  than  1.  A  binomial  series  expansion  can  therefore  be  used 
which  is  formally 


(l+b)l/* 


— i — b2  + 
(2X4) 


(DO).. 

(2X4X6) 


(6) 


The  sin(a)  being  Identified  as  P,  the  expansion  is  applied  to  yield 
separate  terms  for  r  from  eqn  (6).  After  some  algebra, 


X*  t3P  X* 

r  =  R  -  *P  +  -11/2  +  p.|  +  +  0  (7) 

where  Q  represents  terms  with  ratios  of  x“+1/R»  with  n  being  2  or 
above.  Each  term  is  then  examined  to  appraise  its  impact  on  the 
result  of  explikr)  when  substituted.  We  choose  for  example  a  satel¬ 
lite  at  100  km  viewed  from  a  10  cm  optical  system  with  a  wavelength 
at  5xl0-7  meters.  The  maximum  values  for  the  variables  are  applied 
for  each  term. 

1.  R  is  a  constant  and  thus  comes  out  of  the  integral. 

2.  kxP  varies  from  +  to  -  (2ir/5xl0-7)(5xl0-2)  =  6.28x10*. 

3.  Term  3  varies  from  0  to  (2ir/5xl0-7)(6xl0'2)2/(lxl0s).  a 
difference  ofT/10  or  5%  of  a  phase  angle. 

4.  Term  4  varies  from  0  to  Or/5xlO-7)(5xlO-2)3/(lxl010).  a 
difference  of  7.85xlO~B,  which  is  negligible. 

5.  Terms  5  and  6  involve  ratios  of  x“+l/RB  smaller  than  that  of 
term  four.  Their  affect  is  therefore  also  negligible. 

If  we  accept  the  scenario  and  assume  that  the  5%  variation  of 
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term  3  makes  an  Insignificant  contribution  to  error,  then  all  but  the 
first  tvo  terms  can  be  dropped.  Thus 


r  *  R  -  xsin(a)  (8) 

Substituting  the  approximation  for  r  back  into  eqn  (3)  and  bring¬ 
ing  the  constant  term  out  of  the  integral  to  form  a  new  complex  con¬ 
stant  ci  yields 

r/2 

D(x)  =  ci  1 0(a)Z(or)expl-ikxsin(a)|da  (9) 

_  t/2 

If  we  let  q  =  sin(a)/A  ,  then  dq/da  =  cos(a)/A  .  The  limits  thus 
change  to  form  the  expression  (with  a  new  aggregate  constant  C2) 

r1/x 

D(x)  =  c*  I  0(arcsin(qA))Z(arcsin<q\))cos-l(arcsln(qx))expl-2wixq|dq  (10) 

-  lA 

If  we  define  S(arcsin(q\))  =  ciZ(<*)cos~l(<x)  ,  then 

1/x 

D(x)  =  j*0(arcsin(qA))S(arcsin(qA>)exp(-2Tixq|dq  (11) 

-  1/A 

Finally,  If  the  object  and  modified  obliquity  functions  are  transformed 
by  a  "q  transform"  into  functions  in  "q"  space,  such  that  CMq)  = 
0(arcsin(qA))  and  S«(q)  =  S(arcsin(qx))  then 

1/A 

D(x)  =  j*Oq(q)Sq(q)exp(-2Tixq|dq  (12) 

-i/x 
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vhich  of  course  is  simply  the  Fourier  transform  of  the  q  transformed 
object  function  modulated  in  amplitude  by  a  modified  and  q  transformed 
obliquity  function.  Whereas  the  rectangular  object  space  provides 
Fourier  transform  properties  only  in  a  small  object  space  region,  the 
spherical  object  space  results  In  approximate  Fourier  propagation  from 
any  oblect  location.  Also  the  rectangular  space  could  not  accomodate 
energy  below  the  plane  of  the  optical  system  aperture  since  this 
plane  corresponds  to  Infinity  in  that  object  space.  The  spherical 
space  of  course  continues  past  this  plane. 

Field  of  View  versus  Obliquity. 

It  is  quite  obvious  that  if  the  product  0*(q)Sq(q)  is  recovered, 
then  a  simple  division  by  Sq(q)  will  yield  our  goal.  However,  if  the 
original  obliquity  function  Z(<*)  behaves  as  a  result  of  vector  addition 
and/or  polarization  phenomena,  then  It  will  likely  have  something  like  a 
cost®)  dependence,  at  least  near  the  optical  axis.  Since  the  modified 
obliquity  function  S*(q)  is  a  result  of  Z(ar)cos"l(a),  the  assumed  be¬ 
haviour  of  Z<®)  will  cause  a  cancellation  of  the  cosine  functions. 

Thus  for  regions  close  to  the  optical  axis  (for  example,  a  ten  degree 
FOV  (field  of  view)),  the  obliquity  function  is  assumed  to  be  unity. 
Certainly  the  optical  system  is  exposed  to  more  than  this  limited  FOV, 
but  since  the  modulation  of  the  spectrum  only  affects  amplitude,  the 
form  of  the  obliquity  function  outside  the  FOV  of  interest  can  be  ig¬ 
nored. 
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Analysis  of  the  New  Constraint. 


Gerchberg  showed  that  when  the  finite  extent  of  the  object 
function  is  known,  the  superresolution  algorithm  does  Indeed  work  to 
continue  the  spectrum.  He  then  tested  the  algorithm  on  the  same  ob¬ 
ject,  but  set  the  finite  extent  limits  greater  than  the  true  limits. 
Although  the  spectrum  was  still  continued  somewhat,  a  definite  modula¬ 
tion  of  the  continued  portion  was  observed  such  that  the  new 
spectrum  approached  zero  with  higher  positive  and  negative 
frequencies.  Because  the  new  constraint  developed,  le  the  semicir¬ 
cular  space,  is  likely  to  be  larger  than  any  object  function,  this 
same  behaviour  is  to  be  expected.  Obviously  this  artifact  will  have  a 
greater  impact  on  actual  spectra  which  have  significant  energy  out¬ 
side  the  band-limits  of  the  optical  system  (ie  high  frequency 
components). 

Incoherent  versus  Coherent  Illumination. 


The  above  development  applies  to  any  object  function  which 
propagates  as  an  infinite  number  of  spherical  waves.  In  the  case  of 
coherent  illumination,  the  object  function  is  the  complex  amplitude  of 
the  radiation  leaving  the  object.  The  Fourier  transform  is  therefore 
the  direct  superposition  of  those  complex  fields  as  they  arrive  at 
the  optical  system.  In  the  case  of  incoherent  illumination,  the  object 
function  is  the  intensity  of  the  radiation  leaving  the  object.  By  the 
Van  Clttert-Zernlke  Theorem,  the  Fourier  transform  of  the  normalized 
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object  Intensity  distribution  under  incoherent  illumination  is  the  com¬ 
plex  degree  of  spatial  coherence  (6:529).  The  problem  of  a  band- 
limited  spectrum  is  similar  in  both  scenarios.  The  advantage  of  in¬ 
coherent  illumination  is  that  the  object  function  -  the  Intensity  -  is 
always  non-negative.  This  quality  allows  a  further  constraint  for  the 
superresolution  algorithm. 

IV.  SUPERRESOLVED  IMAGING  SIMULATIONS 

The  goal  of  superresolution  is  to  develop  a  technique  to  Improve 
resolution  without  increasing  aperture  size.  Thus,  the  test  of  a  su¬ 
perresolution  technique  shall  be  to  determine  if  it  can  yield  an  image 
which  has  the  resolution  of  a  certain  aperture,  given  the  band-passed 
spectrum  of  a  smaller  aperture.  Specifically,  the  band-passed 
spectrum  from  a  4  centimeter  aperture  shall  be  used  to  attempt  to 
build  a  spectrum  of  a  32  centimeter  effective  aperture.  The 
wavelength  used  is  1.6  mm.  The  incoherent  scenario  is  applied  due  to 
its  more  abundant  natural  occurence. 

Two  one-dimensional  test  objects  are  used.  The  first  is  a 
two-point  target  with  an  offset  center.  The  separation  of  the  two 
points  is  chosen  so  that  the  points  are  unresolved  in  the  image 
produced  by  the  4  cm  aperture,  but  are  certainly  resolved  with  the 
32  cm  aperture.  The  second  target  is  a  wide,  offset  three-bar  tar¬ 
get  somewhat  like  an  Air  Force  resolution  target.  This  type  of  ob¬ 
ject  is  more  representative  of  typical  extended  objects.  For  each 


object  Intensity  distribution  the  Fourier  transform  is  applied  analyti¬ 
cally,  then  the  Hartley  spectrum  is  generated  using  eqn  (1) 
(representing  the  received  complex  degree  of  spatial  coherence).  The 
Hartley  spectrum  Is  then  sampled  over  its  center  four  centimeters  at 
64  equal  spacings.  These  values  are  entered  into  the  center  of  the 
initially  zeroed  612  element  RESULT  array  as  the  initial  algorithm 
spectrum  (512  elements  are  required  for  the  32  cm  aperture).  The  ad¬ 
ditional  constraint  of  a  non-negative  intensity  is  also  applied  in  the 
algorithm. 

The  computer  code  used  is  described  in  Appendix  B.  The  results 
are  shown  in  Figures  3  and  4.  Each  plot  in  the  left  column  is  made 
up  of  the  612  values  of  the  Hartley  spectrum  in  frequency  space. 

The  plots  in  the  right  column  Illustrate  the  image  constructed  from 
the  spectrum  on  its  left.  The  top  set  illustrates  the  diffraction 
limited  spectrum/image  pair  if  an  actual  32  cm  aperture  optical  system 
is  used.  The  center  set  Illustrates  the  spectrum/image  passed  given 
by  the  4  cm  aperture.  The  spectrum  of  this  set  also  represents  the 
starting  point  of  the  algorithm.  The  bottom  set  is  the  result  of  the 
simulation  after  1000  Iterations. 

The  results  of  using  the  two-point  test  object  are  outstanding. 

It  is  apparent  that  with  the  q  space  limits  wider  than  the  object 
boundaries,  the  spectrum  is  modulated  as  predicted.  The  image  of  the 
three-bar  test  object  is  also  remarkably  improved,  although  not  as 
dramatically  as  the  two-point  object  results.  Note  the  actual 


spectrum  has  most  of  its  energy  within  the  band-limited  regime.  Thus 
very  little  error  energy  is  being  subtracted  upon  each  iteration. 

The  explanation  for  a  modulation  when  the  finite  constraints  are 
greater  than  the  actual  object  function  boundaries  is  hypothesized  as 
follows.  The  growth  of  energy  outside  the  band-limited  spectrum  can 
be  considered  to  be  the  result  of  a  transform  of  a  finite, 
reconstructed  image.  The  image  is  finite  due  to  the  truncation  of 
energy  outside  the  semicircular  limits.  Such  a  finite  image  must  have 
a  band-unlimited  spectrum.  Thus  as  the  algorithm  progresses,  the 
energy  "spilled"  outside  the  band  limitations  is  essentially  due  to  the 
diffraction  caused  by  truncating  the  image  function.  The  truncation 
of  error  energy  in  q  image  space  can  be  regarded  as  multiplying  the 
generated  image  by  a  rect  function  the  width  of  the  constraint  in 
that  space.  If  the  constraint  locations  match  the  actual  object 
dimensions,  then  such  a  rect  function  would  have  no  effect.  If 
however,  the  rect  function  is  larger  due  to  a  greater  constraint  such 
as  the  semicircular  constraint,  then  its  affect  can  not  be  ignored. 

The  multiplication  by  the  rect  in  q  space  is  represented  by  a  con¬ 
volution  of  the  image  by  a  sine  function  in  frequency  space.  It  is 
this  convolution  which  allows  energy  to  be  spilled  beyond  the  band- 
limited  spectrum.  Since  the  algorithm  only  truncates  error  energy  in 
q  space,  the  spectrum  resulting  from  the  spilled  energy  must  be  the 
emerging  actual  spectrum.  In  this  view  then,  the  sine  function,  which 
is  convolved  with  the  band-limited  spectrum,  is  responsible  for  deter¬ 
mining  how  much  energy  Is  spilled.  Note  that  the  convolution  is  only 


with  the  band-limited  spectrum.  Thus  the  maximum  spillage  will  occur 
close  to  the  band-limited  spectrum,  as  illustrated  In  the  results. 

As  the  rect  function  representing  the  constraint  in  object  space 
gets  smaller,  the  convolving  sine  function  in  frequency  space  will 
broaden,  allowing  more  energy  to  spill  outside  the  limits  and  increas¬ 
ing  the  growth  of  the  continued  spectrum.  In  the  limit  as  the  rect 
approaches  the  object's  finite  dimensions,  the  spillage  will  actually 
become  the  Fourier  transform  of  the  object  and  error  functions  only 
inside  the  object  function  domain.  Thus  as  the  constraints  get  closer 
to  the  actual  constraints,  the  algorithm  should  converge  quicker. 
Conversely,  with  the  large  rect  function  implied  by  a  semicircular  win¬ 
dow.  the  algorithm  should  converge  slowly  since  the  spillage  is  more 
confined  to  near  the  band-limited  spectrum.  If  this  argument  is  true, 
it  only  implies  a  slowing  of  the  algorithm,  not  necessarily  a  different 
final  result.  Theoretically,  given  enough  time,  the  algorithm  should 
converge  to  the  same  image  regardless  of  how  large  the  constraint 
window  is  (as  long  as  it  is  larger  than  the  actual  window). 

V.  CONCLUSIONS 

The  results  indicate  that  superresolution  is  possible  without  a 
knowledge  of  the  finite  dimensions  of  the  object.  The  requirement  for 
the  technique  is  a  quasi-monochromatic.  band-limited  spectrum.  This 
spectrum  can  be  that  generally  associated  with  any  valid  transform 
relation  between  the  object  space  and  a  frequency  space  which  has 


applicable  constraints.  The  Hartley  transform  relation  provides  a 
working  algorithm  which  is  inherently  faster  and  less  memory  intensive 
than  the  Fourier  relation. 


As  a  final  note,  the  technique  demonstrated  is  not  necessarily 
confined  to  superresolution.  Certain  methods  of  phase  retrieval  are 
also  modelled  after  the  Iterative  approach.  Since  the  examined  con¬ 
straint  is  in  object  space,  any  technique  which  benefits  from  object 
dimension  constraints  will  benefit  from  this  newly  found  constraint. 
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APPENDIX  A 


A  BRIEF  HISTORY  OF  SUPERRESOLUTION 

This  appendix  is  an  account  of  the  development  of  the  field  of 
superresolution.  To  be  complete  in  itself,  some  of  the  development 
from  the  main  thesis  is  duplicated. 

It  must  be  noted  that  many  of  the  attempts  for  superresolution 
are  concerned  with  the  continuation  of  a  known  spectrum,  even  though 
actual  data  only  yields  the  amplitude  of  the  band-limited  spectrum. 
The  recovery  of  phase  to  construct  the  band-limited  spectrum  is  a 
field  in  Itself  and  will  not  be  addressed  in  this  background. 

The  Classical  Limit  of  Resolution. 


Resolution  traditionally  refers  to  the  ability  to  distinguish  be¬ 
tween  closely  spaced  points  in  the  object  space.  Thus  a  resolution 
limit  is  generally  accepted  as  the  minimum  angular  separation  of  two 
object  points,  as  taken  from  the  optical  aperture,  which  yields  an 
image  which  can  be  interpreted  as  that  of  two  points.  We  therefore 
say  that  if  the  angular  separation  is  less  than  this  limit,  the  two 
points  are  unresolvable. 

It  was  Sir  George  Airy  (1801-1892)  who  first  applied  the  Fraun- 
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hofer  diffraction  formula  to  a  distant  point  source  propagating 
through  a  circular  aperture  and  discovered  that  the  field  distribution 
in  the  far  field  is  a  Bessel  function  -  the  intensity  being  later 
called  an  Airy  pattern  (6:419).  In  1896,  Lord  John  William  Strutt 
Rayleigh  generalized  this  interpretation  to  imaging. 

He  envisaged  each  point  on  the  object  as  a  coherent 
source  whose  emitted  wave  was  diffracted  by  the  lens  Into  an 
Airy  pattern.  Each  of  these  in  turn  was  centered  on  the  ideal 
image  point  (on  the  Image  plane)  of  the  corresponding  point 
source.  Thus  the  image  plane  was  covered  with  a  distribution  of 
somewhat  overlapping  and  interfering  Airy  patterns  (6:563). 

With  this  interpretation.  Lord  Rayleigh  rather  arbitrarily  as¬ 
signed  the  limiting  resolution  of  a  system  as  that  angle  between  two 
point  objects  where  the  peak  of  one  Airy  pattern  falls  on  the  first 
minimum  of  the  other. 

In  1966,  G.  Toraldo  Di  Francla  asserted  that  "...from  the  mathe¬ 
matical  standpoint,  the  image  of  two  points,  however  close  to  one 
another,  is  different  from  that  of  one  point  (7:497)."  Upon  this  basis, 
Toraldo  concluded  that  there  is  no  theoretical  limit  of  resolving 
power.  Further,  he  claimed  that  the  only  true  limit  is  a  practical 
one  which  Is  related  to  the  refinements  of  detector  technology.  From 
this  point  on  virtually  every  reference  to  superresolution  acknow¬ 
ledges  this  realization  and  Toraldo's  achievement. 

Approaches  to  Eliminate  the  Practical  Limit. 
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Instead  of  reviewing  the  many  efforts  in  superresolution  in  the 
past  thirty  years  one  by  one,  it  is  more  efficient  to  categorize  them 
along  their  general  approaches  and  examine  the  combined  contribution 
of  each  approach.  The  different  categories  include:  Information 
theory,  analytic  deconvolution.  Iterative  error  reduction  and  decon¬ 
volution  by  matrix  methods. 

Information  Theory.  In  trying  to  create  a  specification  for  the 
practical  limit  of  resolution,  Toraldo  examined  the  number  of  degrees 
of  freedom  which  fully  describe  an  object.  The  value  of  each  of 
these  degrees  is  the  value  of  the  spatial  frequency  distribution  at 
points  given  by  sampling  theory.  He  concludes  that  since  a  true  ob¬ 
ject  spatial  description  must  have  an  infinite  number  of  degrees  of 
freedom,  and  the  image  is  in  fact  limited  in  spatial  frequency  due  to 
the  band-pass  nature  of  the  aperture,  the  complete  object  can  never 
be  recovered.  He  thus  defines  the  limit  of  resolution  as  that 
resolution  corresponding  to  the  maximum  degrees  of  freedom  possible 
for  the  optical  system.  In  fact,  he  claimed  that  there  are  many  dif¬ 
ferent  objects  which  can  have  the  same  narrow  band  of  observed  spa¬ 
tial  frequencies.  He  does  admit  however,  that  any  knowledge  about 
the  object  known  a  priori  must  yield  more  information.  Thus  this  ob¬ 
ject  knowledge  is  the  only  means  of  improving  resolution  beyond 
Toraldo's  limit. 

The  basis  of  Toraldo’s  argument  is  that  the  optical  system 
yields  only  a  portion  of  the  spatial  frequency  content  of  the  object. 
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In  1963  J.  L.  Harris  applied  the  mathematics  of  analytic  functions  to 
Imaging  (3).  He  showed  that  If  an  object  has  a  definite  size,  then  the 
analyticity  of  the  observed  spatial  frequency  function  passed  by  the 
system  can  be  continued  throughout  the  frequency  domain.  Harris' 
conclusions  were  so  important  to  the  future  of  superresolution,  that 
they  are  listed  here  directly  from  his  paper. 


(1)  For  objects  of  finite  angular  dimensions,  knowledge  of 
the  spatial  frequency  spectrum  within  the  passband  of  any  imag¬ 
ing  system  implies  knowledge  of  the  spatial  frequency  spectrum 
over  the  entire  frequency  domain,  and  hence  implies  complete 
knowledge  of  the  object. 

(2)  Two  distinctly  different  objects  of  finite  angular  size 
cannot  have  Identical  images,  so  that  the  ambiguous  image  does 
not  exist  for  realizable  objects. 

(3)  Diffraction,  therefore,  imposes  a  resolution  limit  which 
is  determined  by  the  noise  of  the  system  rather  than  by  some 
absolute  criterion. 

(4)  Object  detail,  absent  In  a  diffraction  image,  may  be  re¬ 
stored  by  a  variety  of  techniques,  with  the  ultimate  limit  in 
such  restoration  Imposed  by  the  system  noise  (3:963). 


To  complete  the  information  theory  approach  to  super¬ 
resolution.  W.  Lukosz  developed  the  techniques  of  exploiting  these 
"limited"  degrees  of  freedom  (8;9).  Lukosz  discovered  that  certain  in¬ 
formation  from  the  image  can  be  enhanced  at  the  expense  of  other  in¬ 
formation.  Although  these  techniques  are  fascinating,  the  above  argu¬ 
ments  by  Harris  show  that  they  are  based  on  a  false  limit  and  there¬ 
fore  may  be  improved. 


Analytic  Deconvolution.  In  accordance  with  Lord  Rayleigh's 
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description  of  an  image,  it  can  be  shovn  that  the  intensity  of  an 
image  point  is  simply  the  square  of  the  convolution  of  the  complex 
object  distribution  and  the  isoplanatic  Impulse  response  of  the  opti¬ 
cal  system  (6:110).  A  technique  to  bring  out  the  object  function  from 
the  convolution  integral  and  place  it  in  terms  of  the  image  function 
is  appropriately  known  as  deconvolution. 

In  1966,  Casper  W.  Barnes,  and  in  1966,  B.  Roy  Prieden 
derived  an  expression  for  the  object  distribution  by  deconvolution 
using  orthogonal  elgenequations  for  objects  of  limited  angular  width 
( 10;  1 1 ).  Noting  that  the  Isoplanatic  impulse  response  for  a  circular 
aperture  is  a  Bessel  function,  they  applied  prior  findings  that  the 
convolution  integral  in  this  case  has  a  complete  orthonormal  set  of 
prolate  spheroidal  eigenfunctions.  Although  the  object  function  was 
brought  out  of  the  integral  successfully,  the  resultant  formula 
remained  a  short  but  Involved  expression  of  Integrals  and  series. 
Prieden  continued  the  treatment  to  the  inclusion  of  noise.  Although 
he  showed  favorable  results  for  certain  types  of  objects,  he  acknow¬ 
ledged  simply  that  the  limit  of  this  type  of  reconstruction  is  noise. 
Similar  conclusions  to  this  effect  were  offered  by  C.  K.  Rushforth  and 
R.  W.  Harris  in  1967  (12).  In  this  paper,  dedicated  to  analyzing  the 
noise  of  the  analytical  deconvolution  method,  it  is  shown  that  as  the 
infinite  series  inside  of  the  object  integral  is  approximated  with  a 
larger  upper  limit  to  reduce  truncation  error,  the  noise  inherent  in 
the  image  becomes  dominant,  yielding  a  trade-off  between  resolution 
and  noise.  The  instigator  of  the  next  method  to  be  discussed  inter- 


preted  these  results  as  nothing  less  than  condemning. 


Rushforth  and  Harris  conclude  that  for  this  type  of  con¬ 
tinuation  (analytic  deconvolution!,  the  increase  in  resolution  for 
a  noisy  diffraction  limited  spectrum,  is  offset  by  the  amplifica¬ 
tion  of  the  error  in  the  reconstructed  function  to  the  extent 
that  achieving  resolutions  much  beyond  the  diffraction  limit  with 
realistically  measured  data,  seems  a  very  dubious  proposition 
(1:71  i). 

Superresolution  by  Iterative  Error  Reduction.  In  1973,  R.  W. 
Gerchberg  developed  an  attractive  superresolution  technique  which 
spawned  a  relatively  great  deal  of  future  work.  An  inverse  Fourier 
transform  of  the  band-limited  spectrum  of  an  object  should  reveal  a 
diffraction  limited  image.  Gerchberg  shows  that  this  diffraction  blurs 
the  boundaries  of  the  image.  If  the  boundaries  of  the  object  are 
known,  then  the  error  associated  with  the  blurred  image  outside  the 
known  boundaries  can  be  reduced.  The  procedure  is  best  explained  in 
his  own  words. 


The  known  portion  of  the  spectrum  is  Fourier  transformed 
to  yield  the  diffraction  limited  (image!  which  is  subsequently 
modified  by  setting  all  of  the  diffraction  limited  (image!  outside 
the  known  extent  of  the  true  object  to  zero.  Thus  modified,  the 
(image!  is  Fourier  transformed  and  the  generated  spectrum  is 
corrected  so  that  that  portion  of  the  frequency  range  in  which 
the  true  spectrum  is  known  is  given  the  correct  values.  The 
procedure  iterates  until  a  correction  criterion  based  on  the  es¬ 
timated  (image!  energy  outside  the  known  extent  of  the  true  ob¬ 
ject  or  the  difference  function  energy  between  the  estimated 
and  true  spectrum  over  the  range  of  the  known  spectrum  is 
satisfied  (1:712). 


Gerchberg  continues  to  show  that  the  error  reduction  per 
iteration  decreases  with  continued  iterations.  He  then  considers 
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noise  as  a  separate  "spectrum"  which  can  be  handled  independently  due 
to  the  linearity  of  the  algorithm  and  concludes  that  since  noise  Is 
not  analytic,  the  algorithm  reduces  the  noise  as  well  as  the  error 
energy  (note:  noise  is  not  be  considered  in  this  thesis,  although  the 
Gerchberg  algorithm  is  fundamental  to  the  techniques  pursued). 

In  1983,  John  G.  Walker  applied  the  Gerchberg  algorithm  to 
coherent  Imaging  and  demonstrated  its  efficacy  by  experiment  (13).  Al¬ 
though  the  results  yielded  an  image  far  better  than  diffraction 
limited,  it  was  not  in  fact  a  true  object  representation.  Walker 
hypothesized  that  discrete  sampling,  imperfect  aperture  edges,  a 
finite  sampling  aperture,  temporal  Instabilities  and  scattering  were 
some  of  the  reasons  for  the  less  than  perfect  result.  Again,  the 
only  a  priori  knowledge  was  the  finite  extent  of  the  object. 

Deconvolution  by  Matrix  Methods.  A  superposition  integral  can 
be  converted  into  matrix  equations  exactly  if  the  matrices  describing 
the  image  and  object  are  infinite  and  the  impulse  response  of  the 
system  is  known.  While  such  infinite  matrices  are  considered,  this 
technique  approaches  analytical  deconvolution  (as  a  special  case). 
Without  considering  the  implications  of  infinite  matrices,  the  super¬ 
position  integral  can  be  reformed  into  the  simple  equation  I  =  DO, 
where  I  is  the  diffraction-limited  image  vector,  0  is  the  object  vec¬ 
tor  and  D  is  a  matrix  of  values  of  the  impulse  response  at  a  position 
defined  by  the  matrix  element  positions  of  I  and  0.  An  interesting 
note  is  that  these  vectors  can  represent  either  intensity  (incoherent 
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propogation)  or  complex  functions  (coherent  propogation).  Solving  for 
the  object  vector  yields 


0  =  D"*I  (A —  1 ) 

The  problem  resides  in  the  fact  that  such  an  inversion 
problem  is  ill  posed  (14:149).  Specifically,  the  calculation  of  D_1  can 
not  by  found  using  traditional  means  due  to  truncation  errors  in  the 
calculation  and  the  presence  of  zeroes  (or  very  small  values)  in  the 
matrix  to  be  inverted.  A  popular  attempt  to  solve  the  problem 
without  inversion  is  by  minimizing  the  reformed  equation  I  -  DO  =  0 
Such  attemps  were  made  by  Crawford,  et  al,  in  1981.  and  again  by 
Stierwalt  in  1985  using  least  squares  methods  (15,16).  Again,  although 
resolution  was  In  fact  enhanced,  warnings  were  given  that  such  a 
technique  is  not  to  be  trusted  in  general,  but  that  a  knowledge  of 
the  object's  extent  will  dramatically  improve  the  results. 

Summary  and  Interpretation  of  Past  Results. 

The  contributions  of  the  four  general  approaches  to  super¬ 
resolution  are  fairly  succinct.  In  the  information  theory  approach, 
the  contribution  is  the  thought  of  trading  off  resolution  of  certain 
characteristics  of  the  object  for  increased  resolution  of  others.  The 
analytical  approach  fails  In  the  presence  of  noise.  The  iterative 
technique  is  successful  only  with  the  knowledge  of  object  boundaries 
and  the  absence  of  radiation  outside  these  boundaries.  It  also 
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presents  an  attractive  insight  on  how  to  reduce  error  In  the  received 
diffracted  distribution.  The  matrix  method  is  very  simple  but  has 
drawbacks  due  to  the  inversion  problem  of  the  diffraction  matrix. 
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APPENDIX  B 


EXPLANATION  OP  COMPUTER  CODE 


The  software  language  used  was  Turbo  Pascal  version  4.0.  The 
computer  used  was  an  IBM  compatible  NEC  multispeed  laptop  running  at 
9.54  MHz,  without  a  math  coprocessor.  The  source  code  Is  given  at 
the  end  of  the  appendix.  Comments  are  supplied  here  Instead  of  in 
the  listing  to  facilitate  following  the  algorithm  within  the  code.  The 
time  required  for  a  512  element  FHT  was  just  over  2  seconds. 


Definitions  and  Termlnolor 


W  The  wavelength  of  incident  radiation. 

width  The  physical  dimension  of  the  array  representing 

the  spectrum  passed  by  the  optical  system. 

numdet  The  number  of  samples  of  the  band-passed  spectrum 

used  In  the  simulation. 

virtual  Those  array  elements  which  correspond  to  samples 

samples  of  frequency  space  outside  the  band-passed 

spectrum. 

totaldet  The  total  number  of  array  elements. 

srllmitlo  The  element  In  q  space  representing  the  left 
boundary  (le  -W/2). 

srlimithi  The  element  in  q  space  representing  the  right 
boundary  (ie  W/2). 

bistart.  The  most  negative  and  positive  band-passed 
bifinish  spectrum  elements,  repectively. 

srstart.  The  most  negative  and  positive  elements  in  the 
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srflnish  array  used  for  superresoiucion,  respectively. 
Description  of  Procedures. 

Define~the_$ystem.  This  procedure  is  where  the  characteristics 
of  the  scenario  are  entered.  For  all  experiments,  the  width  shall  be 
four  cm;  the  wavelength,  1.5  ram;  nuradet.  64  elements;  and  totaldet,  512 
elements.  From  these  initial  parameters,  all  results  can  be  inter¬ 
preted.  The  rather  long  wavelength  is  used  to  provide  scenarios 
which  can  exhibit  diffraction  effects  representative  of  the  Rayleigh 
limit.  Why  this  is  necessary  is  obvious  when  one  considers  the  resul¬ 
tant  sampling  in  q  space.  Obviously  we  desire  enough  data  points  in 
the  result  to  reasonably  depict  the  image  function.  For  example,  a 
sine  function  cannot  be  described  reasonably  by  function  values  at 
only  three  points.  The  more  points  used  to  describe  a  function,  the 
more  the  sampled  function  will  resemble  the  actual  function.  Unfor¬ 
tunately  with  the  FHT  we  cannot  control  the  spacing  of  the  samples  in 
object  space  directly.  The  spacing  between  values  in  q  space  is 
equal  to  the  inverse  of  the  width  between  values  in  frequency  space 
(another  consequence  of  the  discrete  transform  used).  Since  the  size 
of  q  space  where  an  object  is  allowed  is  limited  to  2/W,  it  Is  obvious 
that  only  a  certain  number  of  equally  spaced  values  are  available  to 
represent  the  reconstructed  object.  Another  way  of  viewing  the  ef¬ 
fect  is  that  if  the  image  has  most  of  its  energy  between  two  samples, 
then  it  will  be  very  difficult  to  interpret  the  success  of  the  algo¬ 
rithm,  Of  course  once  a  spectrum  is  achieved,  a  continuous  Hartley 
transform  can  be  performed  to  recover  the  object  function  at  any 
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location  In  q  space.  This  will  not  be  done  In  this  thesis. 


Initialize.  This  procedure  defines  many  of  the  parameters  intro¬ 
duced  in  the  above  list  of  terms. 


Display.  This  procedure  simpiy  outputs  an  axis,  a  center  line, 
the  limits  in  object  space  and  any  function  desired  onto  a  monitor 
display  in  640  X  200  graphics  mode. 


FHT.  This  is  the  Fast  Hartley  Transfoi.n  procedure.  It  is  a 
modification  of  a  source  code  presented  by  O'Neill  (17). 


Three_Bar_Object  and  Two_Point_Object.  These  procedures  calcu 
late  the  spectrum  of  the  three-bar  and  two-point  objects  used  in  the 
experiments. 


Continue  Spectrum.  This  procedure  represents  the  supcrresolu- 
tlon  algorithm.  The  input  to  this  procedure  as  shown  is  the  band- 
limited  spectrum  calculated  earlier. 


Source  Code 

PROGRAM  SUPERRESOLUTION; 

USES  crt,graph3; 

CONST  max_full_slze  =  613;  max_half_size  =  257; 

TYPE 

type_ldc  =  array(-max_haif_slze..max_half_size|  of  real; 
type_2d  =  array(1..2,l..max_full_sizel  of  real; 
dlrection_type  =  (direct, reverse);  dlsplay_type  =  (object. spectrum); 
VAR 

result, htly_BL_spect  :  type_ldc;  perfunc  :  type_2d; 
per_raat  :  array! l..max_full_size|  of  integer; 
totaldet.srllmitlo.srlimithi, bistart, bifinish, smart, srflnish, 


pwr,n2,ct, iteration, numdet  ;  integer; 
incX,trap,rt2, width, wavelength  :  real; 
converged.INITIALIZED  :  boolean; 


PROCEDURE  DEFINE_THE_SYSTEM; 
BEGIN 


width 

:=  4; 

icml 

numdet 

:*  64; 

(must  be  2**pi 

wavelength 

:=  0.15; 

Icml 

totaldet 

:=  512; 

END; 

PROCEDURE  INITIALIZE; 

BEGIN 

rt2  :=  sqrt(2.0); 

srllmitlo  :=  -(trunc(totaldet*width/(numdet‘wavelength))+l); 
srlimithi  :=  -srllmitlo; 

srstart  :=  l-(totaldet  div  2);  srfinish  :=  -srstart*  1: 

blstart  :=  1 -(numdet  div  2);  blfinish  :=  -prstart; 

for  ct  :=  -256  to  256  do 
begin 

htly_JBL_spectlct]:=0.0;  resultlctl:=0.0; 
end; 

IncX  :=  width/numdet; 

END; 

PROCEDURE  DlSPLAY(mat:type_ldc;typ:display__type); 

VAR  k:integer;  inc,mx;real;  ch:char; 

BEGIN 

hires;clearscreen;draw(64, 100,576, 100,1);  inc:=511/(totaldet-l);mx:=0; 

for  k  ;=  0  to  totaldet-1  do 

begin 

If  k=-srstart  then  draw(64+round(k*inc),  130, 64+round(k‘inc),  160,1); 
if  typ=obJect  then  if  ((k+srstart)=srlimltlo)  or 
((k+srstart)=srlimithl)  then 
draw(64+round(k*inc),112,64+round(k*inc),  120,1); 

end; 

if  typ=obJect  then  for  k  :=  srstart  to  srfinish  do 
mat(k);smat(kj*matfk); 

for  k:=*srstart  to  srfinish  do  if  abs(mat(k|)>mx  then  mx:=  abs(matlkl); 
if  mx>0  then  for  k  :*  srstart  to  srfinish  do  matfk|:=80*matfki/mx; 
if  mx>0  then  for  k:=  0  to  totaldet-1  do 

draw(64+round(k*inc),100,64+round(k‘inc),l00- 

round(matlk+srstart|),l); 

END; 

PROCEDURE  FHT(d_ln:type  ldc;  n:integer;  direction:direction_type): 

VAR 

dt  :  type_2d;  ang, omega, mult  :  real; 

c,d,i,j,k,s,t,ta, fa, t_index, modify, s_start,s_end, power.  trig_ind. 
trig_inc,temp_ind,i_temp,tpl, section  :  integer; 

BEGIN 
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IF  INITIALIZED=FALSE  THEN 
BEGIN 

INITIALIZED:=TRUE; 

pwr:=round(ln(n)/ln(2));  ang:=0.0;  omega:=2.0*pl/n;  n2:=  n  div  2; 

for  1  :=  1  to  n  do 

begin 

perfunc(l,i|:ssin(ang);  perfunc(2,i|:=cos(ang);  ang:=ang+omega; 
d:=0;  t_index:=i-l; 
for  c  :=  1  to  pwr  do 
begin 

s  :=  t_index  div  2;  d  :=  d+d+t_index-s-s;  t_index:=s; 
end; 

per_mat(i|:=d+l; 

end; 

END; 

power:=l;  fa:=l ;  ta:=2; 

for  1;»1  to  n2+l  do  dt(fa,per_matli|l:=d_ln(l-ll;for  l:=n2+2  to  n  do 
dt(fa,per_mat(i||:=d_ln(i-n- 1 1; 
for  i  :=  1  to  pwr  do 
begin 

J:=l;  section:-l;  trig_inc  :=  n  div  (power+power); 

while  ]  <=  n  do 

begin 

trig_ind:=l;  s_start:=section*power+l;  s_end:=(section+l) ’power; 
tpl  :*s_start+s_end+ 1 ; 
for  k  1  to  power  do 
begin 

t:=J+power;  if  (s_start=t)  or  (pwr<3)  then  modlfy:=t  else 

modify:=tpl-t; 

dt(tajl:=dt(faj|+dt(fa,tl*perfuncf2,trig_lnd|+ 
dt(fa,modify|*perfuncfl,trig_ind|; 
temp_lnd: =trig_ind+n2 ; 

dt(ta,t):=dt(faJJ+dt(fa,t|,perfuncf2,temp_indi+ 
dt(fa,modlfy|*perfuncfl,temp_ind|; 
trlg_lnd:=trig_ind+trig_inc;  j:=J+ 1 ; 
end; 

j  :=j4-power;  section:=section+2; 
end; 

power:=power+power;  i_temp:=ta;  ta:=fa;  fa:=i_temp; 
end; 

if  direction=direct  then  mult:=1.0*n  else  mult:=1.0; 
for  i  :=  1  to  l+n2  do  resultli-l):=dtlfa,i|/mult; 
for  i  :=  n2+2  to  n  do  resultll-n-1  J:=dt[fa,i|/mult; 

END; 

PROCEDURE  TWO_POINT_OBJECT;  var  al,a2  :  real; 

BEGIN 

al:=0.050;  a2:=0.2; 

for  ct  :=  prstart  to  prflnlsh  do 

begin 

htly_BL_spect(ct|:=cos(2*pi*al*ctMncX/wavelength)+ 

co8(2*pi*a2*ct*lncX/wavelength)-sin(2'pi*al*ct'incX/wavelength)- 
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sin(2*pi'a2*cfincX/wavelength); 

end; 

END; 

PROCEDURE  THREE  BAR  OBJECT;  var  A.B.C  :  real;  amp, phase  :  type_ldc; 
BEGIN 

A  :=0.12;  C:«0.15;  B:=0.23; 
for  ct  :=  prstart  to  prflnish  do 
begin 

if  ct<>0  then 
begin 

amp|ctl:=(sin(pi*A*ct*incX/wavelength)/(pl*A*ct*incX/wavelength))* 
(l+2*cos(2*pl*ct*incX*B/wavelength)); 
end  else  amp(ct]:=3.0; 

phase{ct|:=cos(2*pi*ct*incX*C/wavelength)+ 

sin(2*pl*ct*incX*C/wavelength); 

htly_BL_8pect(ctl:=amp(ctl*phase[ct); 

end; 

END; 

PROCEDURE  CONTINUEJSPECTRUM; 

BEGIN 

INITIALIZED: =FALSE; 

iteration  :=  0;  converged:=false;  , 

while  converged=false  do 

begin 

for  ct  :=  prstart  to  prflnish  do  resultfctl  htly_BL_spectfct); 
FHT(result,totaldet, reverse); 
for  ct  :=  srstart  to  srlimitlo  do  resultlct):=0.0; 
for  ct  :=  srlimithl  to  srflnish  do  result[ct):=0.0; 
for  ct  :=  srlimitlo  to  srlimithl  do  if  result(ctl<0.0  then 
result(ctl:=0.0; 

FHT(result,totaidet. direct);  gotoxy(  1 , 1 ); writeln(lteration); 
iteration:»iteration+ 1 ; 
if  iterations  1000  then  converged;strue; 
end; 

END; 

BEGIN  {  MAIN  PROGRAM  SECTION  I 

define_the_system; 
initialize; 

|two_point_object;  OR}  three_bar_object; 
contlnue_spectrum; 

END. 
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